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A NONLINEAR MIXED EFFECTS DIRECTIONAL MODEL FOR 
THE ESTIMATION OF THE ROTATION AXES OF THE 
HUMAN ANKLE 

By Mohammed Haddou, Louis-Paul Rivest^ and 
Michael Pierrynowski 
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This paper suggests a nonlinear mixed effects model for data 
points in SO (3), the set of 3 x 3 rotation matrices, collected accord- 
ing to a repeated measure design. Each sample individual contributes 
a sequence of rotation matrices giving the relative orientations of the 
right foot with respect to the right lower leg as its ankle moves. The 
random effects are the five angles characterizing the orientation of 
the two rotation axes of a subject's right ankle. The fixed parame- 
ters are the average value of these angles and their variances within 
the population. The algorithms to fit nonlinear mixed effects models 
presented in Pinheiro and Bates (2000) are adapted to the new di- 
rectional model. The estimation of the random effects are of interest 
since they give predictions of the rotation axes of an individual ankle. 
The performance of these algorithms is investigated in a Monte Carlo 
study. The analysis of two data sets is presented. In the biomechani- 
cal literature, there is no consensus on an in vivo method to estimate 
the two rotation axes of the ankle. The new model is promising. The 
estimates obtained from a sample of volunteers are shown to be in 
agreement with the clinically accepted results of Inman (1976), ob- 
tained by manipulating cadavers. The repeated measure directional 
model presented in this paper is developed for a particular appli- 
cation. The approach is, however, general and might be applied to 
other models provided that the random directional effects are clus- 
tered around their mean values. 

1. Introduction. The human ankle joint complex has been modeled as a 
two fixed axis mechanism. It is the primary joint involved in the motion of 
the rearfoot with respect to the lower leg. The characterization of walking 
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TT Deviation Angle TT Inclination Angle 




ST DeviaWon Angle ST Inclination Angle 




Fig. 1. The deviation and inclination angles of the tibiotarsal (TT) and subtalar (ST) 
human ankle rotation axes. 

disorders associated with cerebral palsy, clubfoot or flatfoot deformities uses 
altered external moments (torques) about the two rotation axes of the ankle. 
An accurate and reliable determination of the orientation of these two axes is 
important to successfully evaluate and treat patients with these conditions. 
There is no consensus, in the biomechanical literature, on a noninvasive 
method for estimating the location and orientation of these ankle axes in a 
live individual. 

The two rotation axes of the ankle can be recorded in an RFU coordi- 
nate system where the x-axis points Right, the y-axis goes Forward and the 
z-axis goes Up. Anatomically, plantarflexion-dorsiflexion occurs about the 
tibiotarsal, or tt axis, which is attached to the lower leg, while the subtalar, 
or st axis, attached to the calcaneus, is used for the supination-pronation 
motion of the foot. These two axes are presented in Figure 1. Their orien- 
tations are determined by four anatomical angles (ttinc, ttdev, stinc, stdev) 
giving the inclinations and the deviations of these two axes referenced to 
the RFU coordinate system. 

Using an average generic orientation for each axis results in substantial 
errors [Lewis et al. (2009)] because of an important between subject varia- 
tion in axis locations as characterized in Section 7.1 below. Several in vivo 
estimation methods have been proposed in biomechanical journals [see, for 
instance, van den Bogert, Smith and Nigg (1994) and Lewis et al. (2009)], 
but none was completely successful at estimating the angles in Figure 1. 
The poor numerical results obtained by Lewis, Sommer and Piazza (2006) 
led them to question the validity of the two-axis model for the ankle. 
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The in vivo estimation of the orientation of the two axes of the ankle is 
a statistical problem. The data set for one individual is a sequence of 3 x 3 
rotation matrices giving the relative orientation of the foot of the subject 
with respect to its lower leg as its ankle moves up and down, right and left, as 
much as possible. Rivest, Baillargeon and Pierrynowski (2008) and van den 
Bogert, Smith and Nigg (1994) provide more details about data collection. 

Following van den Bogert, Smith and Nigg (1994), Rivest, Baillargeon and 
Pierrynowski (2008) developed a statistical model for analyzing the rotation 
data collected on a single ankle. Its parameters are the four angles defined 
in Figure 1 and a fifth one for the relative position of the two rotation axes. 
This model fits well; the residual standard deviations reported in Rivest, 
Baillargeon and Pierrynowski (2008) are around one degree. For a subject 
whose ankle has an average range of motion, the estimates are, however, not 
repeatable. Two data sets collected on the same ankle in similar conditions 
can give diff'erent estimated angles. This occurs because the likelihood func- 
tion does not have a clear maximum; it has a plateau and some angles cannot 
be estimated independently of the others. Indeed, Rivest, Baillargeon and 
Pierrynowski (2008) demonstrate that only three parameters can be reliably 
estimated in a subject with an average ankle range of motion. Considering 
the small range of the residual angles, a failure of the two-axis model is an 
unlikely cause for the poor repeatability of the results. 

This paper suggests methods to improve the numerical stability of the 
estimates derived from the ankle model. A penalized likelihood is proposed 
for estimating the parameters. The penalty is obtained by assuming a prior 
multivariate normal distribution for the five angle parameters. When the 
mean and the variance covariance matrix of the prior distribution are not 
known, one is confronted with a nonlinear mixed effects directional model 
whose parameters can be estimated using ankle's data collected on a sample 
of volunteers. Two numerical algorithms to fit this model are proposed. Their 
performances are evaluated in a Monte Carlo experiment; two data sets are 
then analyzed using the new nonlinear mixed effects directional model. 

Nonlinear mixed effects vary with the parametrization of the random 
effects. Thus, the first step of the analysis is to parameterize the model of 
Rivest, Baillargeon and Pierrynowski (2008) in terms of the angles presented 
in Figure 1 and to derive inference techniques using this parametrization. 
These procedures are then generalized to a Bayesian setting obtained by 
multiplying the likelihood by the prior distribution of the model parameters. 
The algorithms of Lindstrom and Bates (1990) [see also Pinheiro and Bates 
(2000)] for fitting nonlinear mixed effects models are then adapted to the 
new directional model. 

Nonlinear mixed effects and Bayesian models with concentrated prior dis- 
tributions could potentially be used in many problems of directional statis- 
tics. These techniques could, for instance, be applied to the spherical regres- 
sion model considered by Kim (1991) and Bingham, Chang and Richards 
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(1992) and to the directional one way ANOVA model of Rancourt, Rivest 
and Asselin (2000) to characterize the between subject variability of the 
mean rotation. 

This new statistical methodology has important applications in biome- 
chanics. By fitting a nonlinear mixed effects directional model to data col- 
lected on a sample of volunteers, one is able to estimate the mean values and 
the between subject variances of the five anatomical angles in the popula- 
tion. These estimates are found to be close to the clinically accepted results 
obtained by Inman (1976) who used direct measurements from cadaveric 
feet. This provides an empirical validation of the two-axis model of van den 
Bogert, Smith and Nigg (1994) for the ankle. The new model also allows the 
analysis of the rotation data collected by Pierrynowski et al. (2003) which 
compared the orientations of the subtalar axis of two groups of individuals 
classified according to their lower extremity injuries. Finally, the penalized 
predictions associated to the mixed effect model are shown to be more sta- 
ble than the estimates obtained by maximizing the standard unpenalized 
likelihood for one subject. 

2. Parameterization of unit vectors and of rotation matrices using anatom- 
iccil angles. Let Ai and B2 be 3 x 1 unit vectors giving the tibiotarsal 
and the subtalar rotation axes in a coordinate system defined according to 
the RFU convention. These unit vectors are first expressed in terms of the 
anatomical angles. Then the Cardan angle decomposition for a 3 x 3 rota- 
tion matrix is briefly reviewed. This section uses the arctan function with 
two arguments, such that arctan(a, b) is the angle whose sine and cosine are 
given by a/Vo^ + and b/y/a'^ + 6^, respectively. 

First consider the tibiotarsal axis, Ai = (^11,^21)^31)^) where de- 
notes a matrix transpose. Formally, the anatomical angles are defined as 
ttinc = ti = — arctan (^31, All) and ttdev = t2 = arctan(A2i) ^ii)- Without 
loss of generality, we assume that the first coordinate of Ai is positive. A 
general expression for unit vectors in the half unit sphere is 



where Dt = a/ 1 + cos^ (t 1 ) tan^ (t2 ) • In a similar manner, one can param- 
eterize the subtalar axis in terms of the anatomical angles stinc = si = 
arctan(i?32, -B22) and stdev = S2 = — arctan(i?i2, i?22) as follows. 



(2.1) 




cos(ti) 
cos(ti) tan(t2) 
— sin(ti) 




(2.2) 
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where Dg = \/l + cos2(si) tan^(s2)- 

The set SO (3) of 3 x 3 rotation matrices is a three-dimensional manifold 
whose properties are reviewed in McCarthy (1990), Chirikjian and Kyatkin 
(2001) and Leon, Masse and Rivest (2006). This paper uses the Cardan angle 
parametrization with the X — Z — Y convention. It expresses an element of 
S0{3) as a function of the Cardan angles a G [— 7r,7r), 7 G [— 7r/2,7r/2) and 
4> G [— vr, vr) as 




— sm7 
cos a cos 7 
sin a cos 7 

where • • • stands for complex trigonometric expressions that are not used 
in the sequel. We also write R = R(a,x) x R(7,z) x R(0,y), where the 
arguments of R(-, •) are the angle and the axis of the rotation respectively. 

The model presented in the next section uses rotation matrices, A{ti,t2) 
and B(si,S2), whose first and second columns are respectively equal to Ai 
and B2. These matrices are given by 

A(ti, t2) = R(ti, y)R[arctan{cos(ti) tan(t2), 1}, -z] 

^ / cos(ti) — cos^(ti) tan(t2) sin(ti)D( 

= — cos(ti)tan(t2) 1 

* y — sin(ti) sin(ti) cos(ti) tan(t2) cos{ti)Dt 

and 

B(si, S2) = R(si, x)R[arctan{cos(si) tan(s2), ^},z] 

^ / 1 — cos(si) tan(s2) 

= — I cos^(si) tan(s2) cos(si) —sin{si)Ds 

Ysin(si) cos(si) tan(s2) sin(si) cos(si)Ds 

where Dt and Ds are defined in (2.1) and (2.2), respectively. 



3. The model for estimating the rotation axes of a single ankle. This 
section expresses the model of Rivest, Baillargeon and Pierrynowski (2008) 
for the estimation of the anatomical angles for a single subject in terms of 
the rotation matrices A(ti,t2) and B(si,S2)- The data set is a sequence of 
time ordered 3x3 rotation matrices {Rj :i = 1, . . . ,n}. The model for Rj 
involves the four anatomical angles, rotation angles {ai:i = l,...,n} and 
{(pi :i = 1, . . . ,n} in [— vr, vr) about the two rotation axes and the angle 70 G 
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(— 7r/2, 7r/2), related to the relative position of the two axes. The predicted 
value for Rj is given by 

(3.1) ^, = A{ti,t2)'R{a„x)-R{-fo,z)-R{(t>i,y)B{si,S2)'^. 

The errors are assumed to have a symmetric Fisher-von Mises matrix 
distribution with density /(E) = exp{K;tr(E)}/cK, E€ 5*0(3), where is 
a normalizing constant; see Mardia and Jupp (2000). If the parameter k is 
assumed to be large so that the error rotations are clustered around the 
identity matrix I3, one has 



(3.2) E: 



where the entries £1, £2, £3 of the skew-symmetric matrix have independent 
M{0, 1/(2k)} distributions; 1/{2k) is called the residual variance. The model 
postulates that Rj = ^jEj, for i = 1, . . . , n. The likelihood is L[ti,t2,si,S2,^0j 
K, {oj}, {(/jj}] = /(^i^Rj). Rivest, Baillargeon and Pierrynowski (2008) 
show that the angles {oi} and {</){} can be profiled out of the likelihood. 
Indeed, 

L[ti,t2, Si, S2,'yo, K, {tti} , {(pi}] < Lp{ti,t2, Si, 32,^0, k) 

(3.3) 

1 r 

= -exp K^{2cos(0f-7o) + l} 

'^^ L i=i 

where 6? = — arcsin(Ai^RiB2) is the Z-Cardan angle of A{ti,t2)~^tli B(si, S2) 
in the X — Z — Y convention; see (2.3). 

Several methods are available to maximize (3.3). However, for the imple- 
mentation of the mixed effects directional model, a closed form expression 
for the score vector for /3 = (ti, t2) ^i, S2, 70)^ is needed. This is derived now. 
Observe that cos(6') = 1 - 2sin2(6'/2), where 0/2 can be assumed to lie in 
the interval [— 7r/2, 7r/2). The log profile likelihood is equal to 

logLp(ti,t2,si,S2,7o,K) = -4Ky^sin^( ' ^ ° j - nlogc^ + 3n«;. 

i=i ^ ^ 

The score for 70 is easily evaluated, viz. 

Q ^ ^ / QZ 7q \ / ~ 7o \ 

— logLp{ti,t2,si,S2,jo,K) = 4K^sin|^ ' ^ ° j ^"^^ 2 " j' 

The score for the four remaining parameters involve the following partial 
derivatives that are evaluated using elementary methods. The property that 
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the partial derivatives of a unit vector and the unit vector itself are orthog- 
onal was used to get the following results: 



d_ 
dh 
d_ 

d 

dsi 
d 

ds2 



Ai 
Ai 
B2 
B2 



sin(ii)tan(t2) . 1 . 

n Ao — t^Aq, 

D2 Dt 

cos(ti){l + tan2(i2)} . 

^ 

sin(gi)tan(s2) p ^ 1 „ 

— ^2 Bi+;d:^3, 



COS si{l + tan2(s2)} 

m 



Bi. 



Since Of = — arcsin(Ai^RjB2), the score for ti is given by 



d_ 

dti 



— logLp(ii,t2,si,S2,7o,«^) 



sm 



i=l 



70 



cos 



70 



Vl-(AiTR,B2)2 5ti 



Ai'''RjB2. 



This can be evaluated using the previous expressions for the partial deriva- 
tives. Repeating this for the other anatomical angles leads to 



d_ 
d(3 



log Lp{f3,K 



sm 



(3.4) 



where 



(3.5) 



■ 70 



sm 



■70 



cos 



_7oY 



d 



2 J dp 



X,; 



4.^ 

i=l 

n 

1=1 

cos{(gf-7o)/2} 

^1 - (AiTR,B2) 

/ - ^^^^m^^ A2-R.B2 - ^ A3-^R.B2 \ 
Df Dt 

cos(ti){l + tan2(t2)} ^ Tp, p. 

Jy2 A2 ±ljh$2 

sin(.0tan(.2i^^Tj,^3^^^^^Tj,,^B3 



■7oj 



cos(si){l +tan2(s2)} 



L>2 



Ai R,Bi 



\/l - (AiTR,B2)2 
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Evaluated at /3 + 5(/3), where (5(/3) is a 5 x 1 vector with entries close to 0, 
the score vector (3.4) is 

-4.f;{sin(t^)x. + ix.X.^5(/3)} 

+ O(||5(/3)f) + O(max[|sin{(0f-7o)/2}|||5(/3)||]). 

One can consider that the last two terms are negligible since the residuals 
{Of — 7o) are small. This is standard in the large k asymptotics used to 
approximate the sampling distributions of estimators in a directional model: 
both (3 — P and the errors ej in (3.2) are assumed to be 0(1/-^/^); see, 
for instance, Rivest and Chang (2006). Equating this to yields a simple 
updating formula. Given its current value /3, the updated value is /3 + (5(/3), 
where 

W = -(^EX.X.^^ 'E{2-n(t^)x.}. 

This calculation can be carried out by regressing the residual vector 
[2sin{(^? — 7o)/2}] on the explanatory variables {Xj}. Theorem 1 of Rivest, 
Baillargeon and Pierrynowski (2008) holds and, as k goes to oo, the max- 
imum likelihood estimator (3 is approximately normally distributed. Once 
the model is fitted, the residual variance 1/(2k) can be estimated using the 
sum of the squared residuals, 

n 

- = -E4sinH(0f-7o)/2}, 

i=l 

where the residual angle Of = — arcsin(Ai'''RjB2) is the Z-Cardan angle of 
A'''RjB in the X — Z — Y convention. Using Grood and Suntay (1983) clin- 
ical interpretation, rotations of angle {Of} occur about a floating axis that 
is orthogonal to both the tt and the st axes. Plots of these angles appear in 
Figure 2 of Rivest, Baillargeon and Pierrynowski (2008); their residual stan- 
dard deviation, y^l/2k, is about one degree. The distribution of the resid- 
uals {[2sin{(0? — 7o)/2}]} is usually approximately normal; the normality 
assumption in (3.2) is not violated for most of the data sets investigated. 
Thus, the proposed model fits well to the ankle data. 

Many individuals have a limited ankle range of motion; the domains for 
angles {oi} and in (3.1) are therefore limited. This makes the estimates 
of the anatomical angles numerically unstable. There can be important dif- 
ferences between the estimates calculated on two data sets collected in suc- 
cession on the same ankle; see Rivest, Baillargeon and Pierrynowski (2008). 
Thus, individual measurements do not allow the estimation of a complete 
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set of anatomical angles. This suggests to borrow strength from other indi- 
viduals and to consider a Bayesian model whose prior distribution penalizes 
extreme parameter values. 

4. A Bayesian ankle model. Assume, for now, that the residual variance 
1/{2k) is known and that the 5x1 vector of anatomical angles (3 is ran- 
dom and has a A/'5(/3o, Sq), where /3o is the average vector of anatomical 
angles within the population and the 5x5 variance covariance matrix XIq 
characterizes the variability of the anatomical angles within the population; 
both are assumed to be known. Elements of (3^ and could be set equal to 
the values of Inman (1976), who studied the variability of these angles. We 
assume that Sq is 0(1/(2k)), thus, there is a fixed 5x5 upper triangular 
matrix Aq such that = ■^o ^(^o ^)/(^'^)' where Aq^ is the inverse of 
A(|^. This section presents an algorithm to derive the mode of the posterior 
distribution of (5 and suggests an approximation for its posterior distribu- 
tion. 

The posterior distribution of (5 is proportional to 



exp 




2 

■ ex.^{-KSSE{(3)}. 



(/3-/3o)^Ao"^Ao(/3-/3o) 



The posterior mode /3 is the value of j3 that minimizes SSE. It can be 
evaluated by adapting the regression algorithm of Section 3 to the Bayesian 
framework. The vector of partial derivatives of SSE with respect to f3 is 
easily derived. 



— = 4 J^sinf^^^ )X, + 2Ao^^oiM - MoJ, 



=1 



where the vector of partial derivatives Xj is defined by (3.5). 

Proceeding as in Section 3, one constructs an algorithm for minimizing 
SSE. The current value /3 is updated to /3 + 5(/9), where 

/ n \ -1 



.1=1 



2sm — Xj 



+ AoTAo(/3-/3o) 
An alternative expression for the updated value is 

/3 + 5(/3)= (f^X^Xi^ + Ao^Ac 
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X 




An approximation to the posterior distribution of the anatomical angles is 
given next. 

Proposition 1. As k^oo, the posterior distribution for 13 satisfies 

^/2^(/3-3)~AA5|o, ^^XiX,^ + Ao^Ao 

where Xj denotes the 5x1 vector of partial derivatives defined by (3.5) and 
evaluated at (3. 

The posterior density of 5 = ^/2k{(3 — (3) is proportional to eyiY>{—iiSSE{(3 + 
(5/\/2k)}. The result is proved by taking a second-order Taylor series ex- 
pansion around 5 = 0. The first order derivatives are null and the variance 
covariance matrix of Proposition 1 is obtained by dropping the o{1/(2k;)} 
terms in the matrix of second-order derivatives. 

We now study some frequentist properties of p. Let /fl*^*^ be the true values 
of the anatomical angles for the individual under consideration. Thus, /3^*-* 
is a realization of the A/'5(/3o)So) prior distribution, such that ||/3^*^ — /9q|| 
is 0(1/\/2k). The difference (3 — is 0(1/\/2k). The leading term of 
this difference consists of a linear combination of individual experimental 
errors and of the penalty associated to the prior distribution. To get a closed 
form expression, one can proceed as in Appendix B of Rivest, Baillargeon 
and Pierrynowski (2008). It suffices to carry out a first-order Taylor series 
expansion of (4.1) in terms of the difference 5{(3) = f3 — and of the 
experimental errors. This yields 

^55i?(/3W+<5(/3)) 

{n n > 

5^ e,xf ) + Xf Xf ^5(/3) + Ao^Ao(5(/3) + /3« - (3,) 
i=l i=l ) 

+ 0{1/k), 

where ej is a AA(0, 1/(2k)) random variable that depends on the error matrix 
Ej, and X^*^ denotes the vector of partial derivatives Xj evaluated at the 
true value /3^*-*, with i?j set equal to ^'j in (3.5). Now (3 corresponds to the 
value of 5{[3) for which (4.2) is null, thus, 

^ = /3W-(^|:xf)xP + AoTAo) ' 
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|E^.xf) + Ao^Ao(/3W -/3o)| +0(1/k) 

X |fj(xP/3W - e.)xf) + Ao^Ao/3o| + 0{1/k). 

This expansion provides an approximation for the prediction error of /3 as 
a predictor of P^^\ The posterior variance covariance matrix of /3'-*^ given 
in Proposition 1 is an estimate of the variance covariance matrix of the 
approximate prediction error. 

5. A mixed model for the simultaneous estimation of several sets of 
anatomical angles. We now have Af subjects and rii, 1 <i < A/, observed 
rotation matrices on each subject. The data set consists of the 3x3 rotation 
matrices {Rjj -.1 = 1,.. .,M;j = 1, . . .,ni}. Let /3j = (tii, t2i, sii, S2i, 7oi)^ be 
the anatomical angles for the ith ankle. As in Section 4, the angles Pi are 
assumed to be random deviates with a five-dimensional normal distribution, 
M^{j3Q, So). The fixed parameters /3o and Tiq are assumed to be unknown. 

In a mixed effects model, the fixed regression parameters and the variance 
components are estimated using a marginal likelihood. To estimate [3q and 
So, we construct such a likelihood using the profile likelihood defined in 
(3.3), rather than the full likelihood for the ankle model. This is acceptable 
since this profile likelihood is also a likelihood, constructed by assuming that 
the angles Of^ — 7oi defined in (3.3) have a centered angular von Mises distri- 
bution with shape parameter 2k. Indeed, since k is large, the distribution of 
2sin{(0?j- — 7oi)/2} is approximately AA{0, 1/(2k;)}. Using this approximation 
in the evaluation of the marginal likelihood Lip gives 



M / N K+5)/2 



(5.1) 



exp|-Kf;4sin2(^^^^)-K||Ao(/3,-/3o)fld/3,. 



Because of a highly nonlinear integrand, the above expression cannot be 
evaluated explicitly. To find the maximum value numerically, we adapt the 
algorithm of Lindstrom and Bates (1990) to the directional ankle model. 

For each i, and for fixed (/3o,Ao), the integrand of (5.1) is maximized 
using the method presented in Section 3. Let (3^ be the maximum value 
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for the ith. sample. Using the first order Taylor series expansion derived in 
Section 3 leads to the following approximation: 

2,in(SL_2i) .2sin(^M_2^ +X./(/3.-ft), 

where the angles Of^ and 7oj and the 5x1 vector of partial derivatives 

Xjj are evaluated at (3.^. Using this approximation, and changing variables 
/3j — /3o = z in the integral, (5.1) becomes 



/ exp<^ -i^yZ^y^ ~ ^ij'^'^ - ^ij'^Pn')'^ ~ '^ll^ozlP > dz. 



n 



(k/vt) 



I + XiAoiAoTXiT|i/2 

X exp{-K(y, - Xi/3o)^(I + X^ Aq ^ Tx,^)-i(y, - Xi/3o)}, 

where Xj is the rij x 5 matrix of partial derivatives for the ith unit, and yj is 
the nj X 1 vector whose j entry is given by Xjj^/3j — 2sin{(^f^- — 7oi)/2}. This 
evaluation of the integral with respect to z follows the argument presented 
in Pinheiro and Bates (2000), Section 7.2. It involves the five-dimensional 
normal density with mean vector (Xj~'"Xj + Aq^ Ao)"^X.j^(yi — Xj/3Q) and 
variance covariance matrix equal to (Xj'''Xj + Aq^ Aq)""*^/ (2k). Thus, Lip is 
approximately equal to the likelihood function for the following linear mixed 
effects model: 



(5.2) 



/ Xi ^ 




/Xi 










X2 








Ixm; 




lo" 











b + e, 



X 



M . 



where y is the (X^n^) x 1 vector of the dependent variable, /3q is the mean 
vector of the anatomical angles in the population, and b is the 5M x 1 vec- 
tor of the individual random effects bj = /3j — /3o that are assumed to be 
independent random vectors with a A/5(0, Sq) distribution. Finally, e is the 
(^nj) X 1 vector of the experimental errors of the directional model con- 
taining independent AA(0, 1/(2k)) random deviates. The second step of the 
algorithm estimates /3o, 1/(2k) and Sq by fitting the linear mixed effects 
model (5.2). This gives updated values for the parameters of the prior dis- 
tribution that are used to get a new set of penalized estimates {/3j}; these 
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new estimates are used to get a new approximation to (5.1) and to update 
the marginal parameter values by fitting (5.2). This two-step algorithm typ- 
ically converges after a few iterations. It is called the PLME algorithm since 
it uses a Penalized least squares and an algorithm for fitting Linear Mixed 
Effects models. Lindstrom and Bates (1990) showed that the first step can 

be bypassed by using fl^ = /3q + as the values around which the lineariza- 
tion of (5.1) is carried out at iteration A; + 1, where is the estimate of 
the random effect for the ith individual at the feth iteration. This one-step 
algorithm is labeled LME. 

This section has considered the one-sample problem, where all the individ- 
uals share the same fixed effect vector /^g. Section 7 considers a two-sample 
model where the direction of the subtalar axis is allowed to vary between 
samples. The directional mixed effects model for this problem is easily con- 
structed. The local linear mixed effects model at step 2 of the Lindstrom 
and Bates (1990) algorithm has a 7 x 1 vector of fixed regression parameters 
featuring 5 entries for the mean angles in the first sample and two param- 
eters for the between group differences of the two subtalar angles. The two 
algorithms proposed in this section can be used to estimate the parameters 
of this enlarged model. 

6. A simulation study. This section reports the results of a Monte Carlo 
experiment to investigate the sampling properties of the estimators of /Bg 
and So obtained by maximizing (5.1) with the two versions of the Lind- 
strom and Bates (1990) algorithm. First, the method used to simulate the 
data is reviewed, then some results will be presented. In the simulations the 
calculations are carried out with angles expressed in radians; for convenience 
the results are presented in degrees. 

Simulations were carried out for the one-sample model only. Values of 
M = 30, 50, 100 and n = 50, 100, 200 were considered. The simulation used 
500 Monte Carlo samples. The following parameter values were used: 



The mean and standard deviations for the first four anatomical angles 
are as given by Inman (1976). The residual standard error 1/\/2k of one 
degree was similar to estimates found in the numerical examples of Rivest, 
Baillargeon and Pierrynowski (2008). 

For each individual, the five anatomical angles were first simulated from a 
A/5(/3o) ^o) and the rotation matrices A{ti,t2) and B(si, S2) were evaluated. 
To construct the predicted values given in (3.1), angles {aij,(pij) obtained 
by fitting the one subject model to some real data were used. The average 
values for {aij,(j)ij) were (38, 14) with standard deviations of (12, 10.5). Thus, 



/3o = (8,-6,42,23,17)^ 



So = diag(7,4,9,ll,ll)2 



(6.1) 



1/(2k) = 1. 
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the motion about the st axis has a smaller range than that about the tt 
axis. The rotation errors Ejj were generated from z = (zi, Z2, z^)'^ three 
independent AA(0, 0.017^) random variables (a standard deviation of 0.017 
radian is 1 when expressed in degrees). Its rotation axis was set to z/||z||, 
while its rotation angle was equal to ||z||. To understand the numerical 
challenges associated to the maximization of the likelihood for the model of 
Section 3, it is convenient to evaluate the vector of partial derivatives Xjj at 
(3 = (3q in an error free model. One gets X^- = (0.01 cos{aij) — 0.99sin(ajj), 
0.99cos(ajj), 0.26 cos((/)ij) + 0.95 sin((/)jj), — 0.80cos((/)jj), 1) . The matrix Xj 
of the vectors of partial derivatives for one subject has a condition number 
larger than 100. This multi-collinearity affects especially ttdev, stdev and 70, 
three rotation angles about different z-axes that are not well differentiated 
when the ankle exhibits a small range of motion. 

The simulations compared the two algorithms proposed by Lindstrom and 
Bates (1990), PLME and LME, as described in Section 5. The R-function 
Ime was used to fit a linear mixed effects model at step 2 of the PLME 
algorithm. This function provides estimates of the sampling variances for the 
fixed effects. The biases of these variance estimators were also investigated 
in the Monte Carlo study. The two algorithms gave almost the same results. 
Only those obtained with PLME are presented. 

Tables 1 and 2 report findings where the fitted model has a diagonal XIq. 
In Table 1, the estimates for ti, t2 and si have small biases, not significantly 
different from 0, and small root mean squared errors. The estimates for S2 
and 7o are less precise. This is caused by the multi-collinearity problem men- 
tioned above. The Ime variance estimates underestimate the true variances; 
this underestimation is more severe for the angles §2, i2 and 70. Increasing n, 
the number of data points by subject reduces this bias. Table 2 is concerned 
with the estimation of the standard deviations. It shows small negative bi- 
ases for all the variances. The parameters and a-.^^ have the largest root 
mean squared errors. Still, Tables 1 and 2 show that the directional mixed 
effects model gives reliable estimates when the random effects are assumed 
to be independent. 

It is likely for the anatomical angles to be correlated so that the true value 
of So niight not be diagonal. Inman (1976) did not consider this question; 
he did not report the correlations between anatomical angles. One could 
investigate this problem by fitting a model with an unstructured Sq with 
15 parameters (that is the variances of the 5 angles and the 10 covariances 
between pairs of angles). Unfortunately, the results obtained with such a 
model are not reliable. In simulations, not reported here, biased estimates of 
the off-diagonal elements of Sq were obtained, especially for the covariances 
involving ^2 1 •S2 and 70 . Apparently, the nonlinear mixed effects model cannot 
distinguish a true correlation between random effects in the population from 
a correlation caused by an ill-conditioned likelihood for the estimation of 



Table 1 

Bias and root mean squared error, in parenthesis, of the estimator of /3q and the relative bias of the Ime variance estimator, in 

parenthesis, when the fitted model assumes that So is diagonal 



n 


M 




Sl 




S2 






7o 






ti 




t2 


50 


30 


-0.09 


(1.70,-2) 


0.67 


(3.10,- 


-16) 


0.55 


(3.27,- 


-14) 


0.10 


(1.56,-12) 


-0.09 


(1.26,-14) 


50 


50 


-0.17 


(1.34,-4) 


0.74 


(2.44,- 


-13) 


0.54 


(2.74,- 


-19) 


0.03 


(1.14,-2) 


-0.06 


(0.99,-12) 


50 


100 


-0.13 


(0.97,-9) 


0.85 


(1.91,- 


-33) 


0.68 


(1.96,- 


-23) 


0.05 


(0.85,-12) 


-0.12 


(0.73,-17) 


100 


30 


-0.09 


(1.58,11) 


0.26 


(2.69,- 


-10) 


0.22 


(2.92,- 


-13) 


0.02 


(1.37,0) 


0.04 


(1.09,-7) 


100 


50 


-0.08 


(1.28,4) 


0.52 


(2.13,- 


-11) 


0.24 


(2.17, 


-4) 


0.06 


(1.06,2) 


-0.03 


(0.79,0) 


100 


100 


-0.07 


(0.90,4) 


0.42 


(1.52,- 


-14) 


0.30 


(1.63,- 


-17) 


-0.09 


(0.79,-8) 


-0.02 


(0.58,-9) 


200 


30 


-0.01 


(1.70,-5) 


0.16 


(2.34, 


-3) 


0.06 


(2.50, 


-9) 


0.06 


(1.38,-7) 


0.05 


(0.86,9) 


200 


50 


-0.09 


(1.32,-5) 


0.39 


(1.87, 


-6) 


0.27 


(1.96, 


-8) 


-0.02 


(1.03,0) 


0.00 


(0.65,13) 


200 


100 


-0.06 


(0.89,4) 


0.21 


(1.34, 


-9) 


0.06 


(1.44,- 


-15) 


-0.01 


(0.75,-5) 


0.05 


(0.53,-13) 



Table 2 

1/2 

Biases and root mean squared errors, m parenthesis, of the estimators of the standard deviations, diag(SQ ), 

when the fitted model assumes that So ts diagonal 



n 


M 




0"t2 


O" SI 


(T. 


82 


(T. 


TO 


50 


30 


-0.05 


(1.10) 


-0.28 


(1.28) 


-0.19 (1.25) 


-0.04 


(2.34) 


-0.46 


(2.50) 


50 


50 


-0.06 


(0.78) 


-0.16 


(0.89) 


-0.12 (0.99) 


-0.02 


(1.74) 


-0.18 


(1.94) 


50 


100 


0.00 


(0.53) 


-0.07 


(0.63) 


-0.05 (0.66) 


0.03 


(1.27) 


-0.10 


(1.26) 


100 


30 


-0.11 


(0.98) 


-0.11 


(0.97) 


-0.07 (1.21) 


-0.14 


(2.09) 


-0.16 


(2.27) 


100 


50 


0.00 


(0.75) 


-0.01 


(0.67) 


-0.06 (0.92) 


0.13 


(1.45) 


-0.02 


(1.66) 


100 


100 


-0.04 


(0.54) 


-0.03 


(0.49) 


-0.05 (0.70) 


-0.02 


(1.12) 


-0.15 


(1.20) 


200 


30 


-0.09 


(0.96) 


-0.06 


(0.76) 


-0.09 (1.14) 


-0.18 


(1.82) 


-0.29 


(1.95) 


200 


50 


-0.05 


(0.75) 


-0.06 


(0.58) 


-0.07 (0.92) 


-0.02 


(1.40) 


-0.12 


(1.51) 


200 


100 


-0.04 


(0.54) 


-0.06 


(0.41) 


(0.68) 


0.03 


(0.97) 


-0.05 


(1.07) 
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Table 3 

The population estimates obtained by fitting the basic model to the 
volunteer data set compared with the values presented by Inman (1976) 







tl 


t2 


Sl 


S2 


70 


Data 




3.32 


-8.05 


38.27 


20.92 


22.42 




s.e. 


0.89 


1.52 


0.89 


2.50 


2.37 




\/ ^Ojj 


5.34 


10.25 


6.46 


14.77 


8.87 


Inman 




8 


-6 


42 


23 


NA 






7 


4 


9 


11 


NA 



the random effects. Additional investigations of this problem could consider 
models where has a small number of nonnull covariances. 

7. Numerical examples. This section presents the analysis of two data 
sets collected in the Human Movement Laboratory of the School of Rehabil- 
itation Science at McMaster University, using an OptoTrak camera system 
at a frequency of 50 Hz; see Rivest, Baillargeon and Pierrynowski (2008) for 
a detailed description of the data collection protocol. The successive rota- 
tion matrices in a data set are not independent measurements; the residual 
autocorrelation when fitting the model of Section 3 on the data collected on 
one subject is larger than 0.80. In order to satisfy, at least approximately, 
the assumption of independence underlying the construction of the penal- 
ized likelihood in Section 4, the model was fitted to a subsample of the data 
obtained by keeping one observation out of 30. A sampling frequency of 1.67 
Hz yielded smaller residual autocorrelations, in the range —0.3 to 0.5; the 
assumption of independence was approximately satisfied. Subsampling the 
data was also used by van den Bogert, Smith and Nigg (1994) to get stable 
estimates of the anatomical angles. 

7.1. An empirical validation of the estimates of Inman (1976). The mean 
values and the population standard deviations of the four anatomical an- 
gles characterizing the direction of the two rotation axes of cadaver ankles 
were presented by Inman (1976), by manipulating unloaded cadavers feet. 
Inman's estimates are given in Table 3. The model of Section 5 provides an 
in vivo method for estimating the same angles. In this section, we use right 
foot data collected on M = 65 volunteers with sample size n = 50 rotation 
matrices. The estimates obtained by fitting the model of Section 5 to this 
data set are also presented in Table 3. This model has an estimated residual 
standard error, l/\/2^ of 0.021 radians, or 1.21 degrees. 

Although two of four of Inman's values are outside of the 95% confidence 
interval, the agreement between the two sets of estimates is reasonably good 
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Fig. 2. Boxplots for two sets of estimates for stinc expressed in degrees. 



(see Table 3). The standard deviations are remarkably close, except possibly 
for the angle t2 = ttdev. Overall, Inman's estimates and the one derived 
using the directional mixed effects model with unloaded ankle motion data 
are similar. 

To investigate the Bayesian model of Section 4, two sets of estimates of 
the anatomical angles of the right ankle of the 65 experimental subjects 
were calculated. The first set used n = 50 observations per subject and the 
Bayesian penalty using (6.1) as the parameters for the prior distribution. 
The second set was obtained by fitting the unpenalized model of Rivest, 
Baillargeon and Pierrynowski (2008) to the complete data sets of n = 1500 
frames per subject. For the two sets of estimates, the mean values for the 
five angles were similar to the estimate 0q in Table 3. There were important 
differences in the between subject standard deviations: for the first set they 
were (5.15,6.62,6.61,12.84,8.25), respectively, for (ii, t2, -si, S2, 7o)) while for 
the second it was (14.72,18.70,15.94,36.56,28.54). Thus, the n = 1500 esti- 
mates are much more variable than their penalized alternatives. The added 
variability is caused by the numerical problems in maximizing the likelihood 
of the unpenalized model. Similar numerical problems were encountered by 
van den Bogert, Smith and Nigg (1994) when they fitted the ankle model 
using an ad hoc loss function. They proposed setting S2 = to get stable es- 
timates. The penalty of the Bayesian model allows to get meaningful values 
for the five angles. 

Several studies have found that si = stinc is the only angle that can be 
reliably estimated using either the model of Rivest, Baillargeon and Pier- 
rynowski (2008) or the approach of van den Bogert, Smith and Nigg (1994). 
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Figure 2 gives the boxplots of the two sets of estimates for si; an outher with 
a negative estimate for si has been removed. The penalty protects against 
values of si = stinc larger than 70 degrees that are not possible anatomically. 

7.2. Data analysis for the location of lower extremity injuries. Pierrynowski 
et al. (2003) collected ankle motion data from 31 participants who had expe- 
rienced either knee (ni = 15) or foot (n2 = 16) injuries during weight-bearing 
activities. The experiment's goal was to determine whether there were sig- 
nificant differences in the orientation of the subtalar axis between the two 
groups. This section uses the data collected on the right ankle; the individual 
sample size is nj = 50 for each subject. 

Table 4 presents the estimates obtained by fitting three models to this 
data; all the models have five independent random effects, one for each 
component of f3. The first one has 7 fixed parameters; in addition to f3, 
it features parameters dsi and ds2 for the foot- knee differences of the two 
subtalar angles. Model 1 postulates that the mean stinc and stdev values are 
(si, S2) and {si + dsi,S2 + ds2) in the knee and the foot group respectively; 
the variances of these two angles are the same for both groups. Models 
2 and 3 are derived from model 1 by setting ds2 = and dsi = ds2 = 0, 
respectively. Under model 3, the mean values of the five angles of the ankle 
model are the same in both groups. The LME and the PLME algorithms gave 
slightly different numerical results; the PLME estimates are presented as a 
larger maximum for the likelihood Lip was obtained with this algorithm. 
The estimated residual standard error \/^/2k was 0.024 radians, or 1.36 
degrees for the three models. 

In model 1, the z-statistic for testing the null hypothesis of no stdev effect 



is Zobs = —2.89/2.56 = — 1.13 for a p- value of 0.26. This high p-value suggests 

Table 4 

The estimates obtained by fitting three models to the lower extremity injury data 

Model ti t2 si dsi S2 ds2 70 

1 /3j -4.14 7.09 45.93 -7.41 6.69 -2.89 -7.33 

1 s.e. 1.49 1.99 1.05 1.42 3.07 2.56 3.22 

1 y^Soji 5.86 6.15 3.68 NA 4.68 NA <10"^ 

2 4 -5.46 7.62 45.78 -7.45 3.44 NA -8.89 
2 s.e. 1.41 1.78 0.99 1.39 2.25 NA 2.54 

2 \/Sojj 5.69 5.40 3.58 NA 3.61 NA <10"^ 

3 4 -4.71 5.94 43.01 NA -0.04 NA -9.75 
3 s.e. 1.49 1.97 0.96 NA 2.69 NA 2.67 

3 a/Soj-j 5.89 7.63 4.88 NA 6.35 NA 4.76 
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setting ds2 = 0; this leads to model 2 where the test for a si between group 
difference has a p-value smaller than 10~^; this is highly significant, even 
when accounting for the underestimation of the standard errors highlighted 
in Table 1. So model 2 provides the best fit, in agreement with the findings 
of Pierrynowski et al. (2003) who also noted the significant difference in 
stinc between the two groups. The nonlinear mixed effects model allows to 
test for an S2 = stdev difference; this was not possible with the unpenalized 
estimates since they were numerically unstable. 

8. A reduced model. In his seminal work, Inman dealt with four angles, 
ttinc, ttdev, stinc and stdev. The fifth angle 70 of (3.1) does not play any 
role in his investigations. This section suggests a reduced model featuring 4 
anatomical angles instead of 5. It investigates whether this model leads to 
better estimates of the anatomical angles. 

In the reference position the leg and the foot reference frames have the 
same orientation, thus, the predicted value 'J' = I is possible. Therefore, for 
some angles a and (p, one has 

1 = A(ti,t2)R(a,x)R(7o,z)R(0,y)B(si,S2)^. 

This equation implies that 70 is the Z-Cardan angle of A{ti,t2)~^B{si, S2) 
in the X — Z — Y convention. Thus, 70 can be expressed in terms of A{ti , ^2) 
and B(si, S2) as 

(8.1) 7o(Ai,B2) = -arcsin(A7B2). 

Fitting a reduced model where 70 depends of (51,52,^1,^2) through (8.1) is 
easily carried out in the framework of Section 3. All the previous derivations 
hold with the reduced model provided that the matrix X is redefined in such 
a way that its ith row is given by the 4x1 vector: 

Xi = -cos{(0,^-7o)/2} 

/_( sinit,)tMt,) j^ ij^y ( ^—]b2\ 

\ Df ^ ' Dt j \^Y^l-(AiTR;B2)2 J 

cos(ti){l+tan2(f2)} A t( Rj I 

A T / Rj I \ / sin(si)tan(s2) p i 1 "R \ 

^1 V vi-(Ai^R.B2)^ ~^jV ^i + iTlBsj 

_ cos(.i){l+tan^(.2)} ^ T/^ R^ 

\ ^? \^vi-(Ai^R»B2)2 ^^-yoy / 

The reduced model has been applied to the data on the 65 volunteers pre- 
sented in Section 7.1. It yielded poor results; the within subject variability 
was larger than that with the five-parameter model of Section 4. The average 
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values failed to reproduce Inman's results. To investigate this failure, we cal- 
culated the differences 70 + arcsin(A^B2) obtained with the five-parameter 
model for the 65 subjects of Section 7.1. The average difference is —2.30 
degree (s.e. = 0.34); only 7 of the 65 differences have positive values. 

The assumption underlying the reduced model, that the leg and the foot 
are aligned in the reference position, is not met. The reference position is 
measured when the experimental subject is standing, so that its ankle is 
loaded. However , the data is collected on an unloaded ankle moving freely 
when the subject is sitting. The nonnull value of 70 + arcsin(A^B2) might be 
explained by a slight change of the relative orientation of the two reference 
frames when the ankle goes from a loaded to an unloaded position. Thus, the 
reduced model is not suitable for the data analyzed in this work. It should, 
however, be considered when investigating the rotation axes of a loaded an- 
kle. The value of 70 -|- arcsin(A7B2) may quantify rearfoot flexibility which 
is of interest to foot care professionals [Mansour et al. (2007)]. 

9. Discussion. This work presented a solution to the estimation of the 
directions of the two rotation axes of the ankle. The key element is the es- 
timation criterion given by a penalized likelihood. This penalized likelihood 
is associated to a Bayesian model for the ankle joint and to a nonlinear 
mixed effects directional model that allows estimation of the between ankle 
variability of the rotation axes within a population. Simulations have shown 
that the population means and the population variances can be estimated in 
a reliable way. When used on a data set collected on a sample of volunteers, 
the nonlinear mixed effects directional model produced mean and variance 
estimates that were similar to those presented by Inman (1976). The good 
match with Inman's clinically accepted findings (see Table 3) provides em- 
pirical evidence that a two-axis (revolute) mechanistic model of the ankle 
[see equation (3.1)] is indeed appropriate for the ankle. Section 7 shows that 
the proposed nonlinear mixed effects directional model can be extended to 
compare the ankle's axes in several populations. The Bayesian model of Sec- 
tion 4 might very well solve the problem of estimating in vivo the location 
of the ankle's rotation axes. 

Future work of this model includes a detailed investigation of the within 
subject stability of the Bayesian estimates of Section 4, using both right and 
left foot data. The estimation of the translation parameters of the van den 
Bogert, Smith and Nigg (1994) ankle's model will also be investigated and 
the application of the reduced model to data collected on a loaded ankle will 
also be considered. 
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